home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 46 / Amiga Format CD46 (1999-10-20)(Future Publishing)(GB)[!][issue 1999-12].iso / -in_the_mag- / reader_requests / scilab / demos / wheel2 / maple / euler.map < prev    next >
Text File  |  1999-09-16  |  5KB  |  182 lines
  1. #----------------------------------------------------------------------------
  2. #     GENERAL PART
  3. #----------------------------------------------------------------------------
  4. #----------------------------------------------------------------------------
  5. # Notation for Tex output 
  6. # and TeX output routines
  7. # this routines generates proper derivative notations up to second derivative
  8. # for a set of variable. The TeX name for the variables is also given
  9. # Warning this function will detroy previous values of texsub 
  10. # Ex : addnotations( { x = `x`  , th = `\\theta`, y =`y` } );
  11. #-----------------------------------------------------------------------------
  12.  
  13. addnotations:=proc(varlist)
  14.     local x,nots:
  15.     nots:={};
  16.     for x in varlist do
  17.         nots:={op(nots),x, 
  18.             cat(lhs(x),`d`)=cat(`\\dot{`,rhs(x),`}`),
  19.             cat(lhs(x),`dd`)=cat(`\\ddot{`,rhs(x),`}`)}
  20.     od;
  21.     texsub:=nots
  22. end:
  23.  
  24. sortieinit:=proc(filename)
  25.     writeto(filename);
  26.     writeto(terminal);
  27. end:
  28.  
  29. sorties:=proc(filename,comment,exp)
  30.     appendto(filename);
  31.     lprint(comment);
  32.     lprint(` `);
  33.     mtex(exp);
  34.     writeto(terminal);
  35. end:
  36.  
  37. sortiesI:=proc(filename,comment)
  38.     appendto(filename);
  39.     lprint(comment);
  40.     writeto(terminal);
  41. end:
  42.  
  43. sortiesM:=proc(filename,comment,expr)
  44.     appendto(filename);
  45.     lprint(comment);
  46.     lprint(` `);
  47.     map(x -> mtex(x),expr);
  48.     writeto(terminal);
  49. end:
  50.  
  51. readlib(tex):
  52. texwid:=120;
  53. sortieinit(`systeme.tex`);
  54.  
  55. #---------------------------------------------------------------
  56. # functions for computing euler equations 
  57. # L : the Lagrangian,
  58. # q,qd,qdd are the state vector and its derivatives
  59. #---------------------------------------------------------------
  60. with(`linalg`):
  61.  
  62. euler_equations:=proc(L,q,qd,qdd)
  63.      local k,m:
  64.     m:=nops(q);
  65.     v:=matrix(m,1,0);
  66.     for i to m  do  
  67.            v[i,1]:=LL(q[i])=simplify(time_diff(diff(L,qd[i]),q,qd,qdd)
  68.                 -diff(L,q[i]));
  69.     od;
  70.     eval(v):
  71.     end:
  72.  
  73. #---------------------------------------------------------------
  74. # Time derivative computation of an expression 
  75. # depending on q,qd,qdd 
  76. # used to compute Euler equations 
  77. #---------------------------------------------------------------
  78.  
  79.     ttvar:=proc(xx)
  80.             if type(xx,`indexed`) 
  81.         then cat(op(0,xx),`d`)[op(xx)]
  82.         else cat(xx,`d`) fi
  83.     end:
  84.  
  85.  
  86. time_diff:=proc(phi,q,qd,qdd)
  87.     local phi_copy,k,ttvar:
  88.     # subtitution to specify that q,qd ,qdd depends on time 
  89.     phi_copy:=phi:
  90.     phi_copy:=subs(map( xx-> xx=xx(t),[op(q),op(qd)]),phi_copy):
  91.     diff_phi:=diff(phi_copy,t):
  92.     # subtitution to come back to our variables 
  93.     diff_phi:=subs(map(xx->diff(xx(t),t)=ttvar(xx),[op(q),op(qd)]),
  94.             diff_phi):
  95.     diff_phi:=subs(map(xx->xx(t)=xx,[op(q),op(qd),op(qdd)]),diff_phi):
  96. end:
  97.  
  98.  
  99. #-----------------------------------------------------
  100. # Rewritting the Euler equations to have a canonical form 
  101. #           ..         .
  102. # El= ME(q)  q + CE(q) q^2 + RE(q)
  103. # Computation of ME,CE,RE 
  104. # CEuler returns a list [ME,CE,RE];
  105. #-----------------------------------------------------
  106.  
  107. CEuler:=proc(E,q,qd,qdd)
  108.     local Me,Ce,Re:
  109.     Me:=MME(E,q,qd,qdd):
  110.     Ce:=CCE(E,q,qd,qdd):
  111.     Re:=RRE(E,Me,Ce,q,qd,qdd):
  112.     [eval(Me),eval(Ce),eval(Re)]:
  113.     end:
  114.  
  115.  
  116. #-----------------------------------------------------
  117. # Rewritting the Euler equations to have a canonical form 
  118. #           ..      .
  119. # El= ME(q)  q + RE(q,q)
  120. # Computation of ME,RE 
  121. # CEulerP returns a list [ME,RE];
  122. #-----------------------------------------------------
  123.  
  124. CEulerP:=proc(E,q,qd,qdd)
  125.     local Me,Ce,Re:
  126.     Me:=MME(E,q,qd,qdd):
  127.     Ce:=matrix(nops(q),nops(q),0):
  128.     Re:=RRE(E,Me,Ce,q,qd,qdd):
  129.     [eval(Me),eval(Re)]:
  130.     end:
  131.  
  132. MME:=proc(E,q,qd,qdd)
  133.     local E1:
  134.     E1:=eval(E):
  135.     genmatrix([seq(E1[i,1],i=1..nops(qdd))],qdd):
  136.     end:
  137.  
  138. #-----------------------------------------------------
  139. # extract the CE(q) matrix  El= ME(q)  qdd + CE(q) qd^2 + RE(q)
  140. #-----------------------------------------------------
  141.     ttvarp:=proc(xx)
  142.             if type(xx,`indexed`) 
  143.         then cat(op(0,xx),`2`)[op(xx)]
  144.         else cat(xx,`2`) fi
  145.     end:
  146.  
  147. CCE:=proc(E,q,qd,qdd)
  148.      local E_copy,q2d:
  149.     E_copy:=eval(E):
  150.     q2d:= map(x-> ttvarp(x),qd): 
  151.     E_copy:=subs( map(x-> x**2=ttvarp(x),qd),eval(E_copy));
  152.     genmatrix([seq(E_copy[i,1],i=1..nops(qd))],q2d);
  153.     end:
  154.  
  155. #-----------------------------------------------------
  156. # extract the RE(q) matrix  El= ME(q)  qdd + CE(q) qd^2 + RE(q)
  157. #-----------------------------------------------------
  158.  
  159. RRE:=proc(E,ME,CE,q,qd,qdd)
  160.     local MM:
  161.     MM:=add(    E,
  162.         add(multiply(ME,matrix(nops(q),1,qdd)),
  163.              multiply(CE,matrix(nops(q),1,map(x->x**2,qd)))),
  164.         1,-1);
  165.     MM:=map((x)-> simplify(x),eval(MM)):
  166.     end:
  167.  
  168. #
  169. #-----------------------------------------------------
  170. # FORTRAN GENERATION 
  171. #-----------------------------------------------------
  172. read ``.libname.`/macrofort.m`;
  173.  
  174. Gener:=proc(filename,fortranlist)
  175.     init_genfor();
  176.     optimized:=true:
  177.     writeto(filename):
  178.     genfor(flist):
  179.     writeto(terminal):
  180.     end:
  181.  
  182.